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ABSTRACT 

The growing field of exoplanetary atmospheric modelling has seen little work on stan- 
dardised benchmark tests for its models, limiting understanding of the dependence 
of results on specific models and conditions. With spatially resolved observations as 
yet difficult to obtain, such a test is invaluable. Although an intercomparison test 
for models of tidally locked gas giant planets has previously been suggested and car- 
ried out, the data provided were limited in terms of comparability. Here, the shallow 
PUMA model is subjected to such a test, and detailed statistics produced to facilitate 
comparison, with both time means and the associated standard deviations displayed, 
removing the time dependence and providing a measure of the variability. Model runs 
have been analysed to determine the variability between resolutions, and the effect of 
resolution on the energy spectra studied. Superrotation is a robust and reproducible 
feature at all resolutions. 
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1 INTRODUCTION 

Over the past several years, work has been carried out in 
the field of three-dimensional numerical simulation of exo- 
planetary atmospheres using large-scale global climate mod- 



els (GCMs), as reviewed in, e.g., Showman, Menou, & Cho 



(2008). In such modelling, it is important to establish clear 



model benchmarks, such as that set up for Earth by " Held fc] 
Suarez ( ^1994 ). Little benchmark work has been carried out 
for exoplanetary studies, however, with the exception of that 



initially described by Menou & Rauscher (2009), who laid 
out an intercomparison test and investigated the response of 
the spectral IGCM dynamical core to it, and further studied 



by Heng, Menou Phillipps ( (2011a ), who investigated the 
response of both spectral and gridpoint cores available for 
implementation in the FMS model to both this and additional 
tests. This study follows the lead provided by these two pa- 
pers in carrying out a clear and reproducible intercompar- 
ison test, and adds to it further statistics and information, 
including online data, to facilitate model comparisons. 

Although many modelling studies have been carried out 
on the atmospheres of 'hot Jupiters', gas giant planets less 
than about 0.1 AU from their host star (e.g. | Johnson|2009 ) , 
few direct comparisons of model responses to these condi- 
tions have been made. It is imperative that model-dependent 
responses be identified and understood. Without such inter- 
comparison studies, it cannot be determined which elements 
of a simulation may correspond to the planet under study, 



and which are simply artefacts of a specific model, or a re- 
sult of the high sensitivity of these complex, non-linear sys- 
tems to the precise input parameters and initial conditions. 



Thrastarson & Cho (2010) have studied the effects of initial 



flow on the final state reached by their model, at atmospheric 
depths from 1 bar down to 100 bar, and found the final state 
to be highly dependent on the initial conditions chosen. In 



contrast, Liu & Showman (2012) carried out a similar study 



down to 200 bar, and found almost no dependence on the ini- 
tial conditions, a result that Liu & Showman (2012|) suggest 



may be due in part to the different vertical profiles of the 
thermal forcing and the restoration timescale between the 
two studies. With little observational information available, 
it becomes vitally important to have a reference simulation 
against which all models can be compared. Whereas models 
of Solar System planets may be compared to observations 
of the planet in question to determine their effectiveness, 
the modelling of extra-solar planets must proceed from first 
principles, though in some cases phase curves may be utilised 
to gain information on the conditions at the photosphere. 

A wide variety of models are utilised in simulating 
hot Jupiter atmospheres, from 'shallow- water' models such 
as that 



used 



Showman & Polvani (2011) to three- 



dimensional GCMs of varying complexity and underlying 
assumptions. Atmospheric depths studied range from 1 bar 
('shallow') to 100s of bars ('deep'), and model atmospheric 
equations range from the one-layer shallow-water equations 
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through the shallow, hydrostatic primitive equations to the 
full three-dimensional Navier- Stokes equati ons (e.g. Dobbs- 
Dixon Lin||2008 ). The SPARc/MiTgcm of Showman et al 



(2009) couples a correlated- /c radiative transfer implementa- 



tion to their dynamical core, while other models may utilise 



dual-band 'grey' radiative transfer (e.g. Heng, Frierson & 



|Phillipps 2011b), or omit radiative transfer and utilise in- 
stead Newtonian relaxation towards a predetermined tem- 
perature state. 

Large variability between models shows that some as- 
pects of the situation have not as yet been positively de- 
termined using such models in their present form. Since 
the input parameters of even simple GCMs are poorly con- 
strained by the available data on the planets to which they 
are applied, the choice of parameter adds a further degree 
of uncertainty to the already variable results. This again il- 
lustrates the importance of a single, fully-specified test case 
from which differing model-dependent responses may be de- 
termined. 

This paper builds on work sugg ested by |Menou fc] 



Rauscher (2009) and further studied in |Heng et al. (2011a), 



and adds the PUMA (Portable University Model of the Atmo- 



sphere) model ( Fraedrich et al.||2005 ) to those to which the 
test has been applied. Since the precise state of the model 
atmosphere is highly time-dependent, long-term statistics 
are produced to allow a quantitative comparison for future 
tests. Time means of the output fields are produced to gain 
an understanding of the overall characteristics of the mod- 
elled atmosphere, and the associated standard deviations to 
gain a quantitative understanding of the degree of variability 
in each field. These statistical analyses of the model output 
fields will allow future work to be compared on a more solid 
footing. 

Section [2] outlines the model used and the description 
of the experiment, including a full list of model parameters. 
Section [3] displays the results of carrying out the intercom- 
par ison test with PUMA, and provides time mean and vari- 
ability statistics. The data files required to produce these 
plots are provided online. Section [4] provides a discussion 
of these results and comparison to the previous results of 
Menou Rauscher] ( [2009t and |Heng et al.| ( |2011al ). Sectionjs] 



covers the conclusions drawn, recommending the use of time 
mean and standard deviation statistics and a minimum sam- 
pling frequency, noting the degree of correspondence to the 
results of the previous authors, and suggesting additional 
plots to capture further aspects of the simulation not seen 
when studying only wind and temperature fields. In com- 
mon with previous studies, and as expected from the work 
of Showman & Polvani (2011), a strong equatorial superro- 



t at ion is found, resulting in an offset temperature hotspot. 



2 A MODEL FOR GAS GIANT EXOPLANETS 

PUMA is a simple global circulation model of the atmosphere, 
developed in its current form at the University of Hamburg 
( Fraedrich et al.|2005 ) and descended from the spectral code 
of Hoskins & Simmons (1975). Though designed for use on 



Earth, the basic equations of the model are applicable to any 
system about which similar assumptions can be made, mod- 
elling an atmosphere in hydrostatic equilibrium that can be 
assumed to be a 'thin shell' with respect to the radius of the 



(spherical) planet. These assumptions, notably the hydro- 
static balance and shallow atmosphere approximations, re- 
sult in a set of equations typically referred to as the primitive 
equations of meteorology. Its long heritage means that the 
model's dynamical core is well known and has been exten- 
sively tested, rendering it particularly useful for benchmark 
tests. 

The model runs in spectral space in the horizontal using 
a triangular truncation, resolution being specified by the 
truncation coefficient, such as T21, which equates to a grid 
of approximately 64 longitudes and 32 latitudes. In general, 
if the truncation number is Tn, the number of latitudes is 
given by A^iat = (3Tn + l)/2. Particular resolutions, such as 
T42 (128x64), T63 (192x96), or T85 (256x128), are widely 
used as standard resolutions owing to the resulting grid sizes 
being entirely or almost entirely powers of two, which allows 
for greater efficiency in the fast Fourier transform routines. 

A finite- difference method is used in the vertical, with 



vertical levels defined using the sigma-coordinate (Phillips 
[1957 ): 



p(A,0,2:,t) 
Ps(A,(/),t) 



(1) 



where p is the pressure at any point in the atmosphere, and 
Ps is the 'surface' pressure at the same A, 0, t, with A, 0, 
z, and t being longitude, latitude, height above the model 
base, and time, respectively, a is then always 1 at the base 
of the model atmosphere, decreasing upwards. Model levels 
(7 = /n, on which the variables are calculated, take values in 
the range < /n < 1, with boundary conditions imposed at 
(7 = and cr = 1. In this gas giant case, the 'surface' is flat, 
determined by a reference geopotential, and the boundary 
conditions equate to a rigid surface at the top and bottom 
of the atmosphere, at which the vertical velocity is set to 
zero. Any level distribution may be input; however, for the 
purposes of this experiment, only linear spacing in sigma 
was used. 

Choosing a horizontal and vertical resolution requires 
a compromise to be made between the degrees of freedom 
of the model and the computing power required to run it 
within a feasible time-scale. Too low a resolution may result 
in misleading results as features vital to the flow are not 
represented. The resolution dependence of this experiment 
is discussed further in Section [4] The particular resolutions 
used in this paper are chosen to correspond with those of 
Menou k Rauscher (2009) and Heng et al^(|2011a). 



PUMA has been modified for use with gas giant plan- 
ets, permitting all planetary parameters to be changed, and 
also to permit the usage of three-dimensional, customised 
temperature restoration fields. Equations [2] to |5] were added 
as a selectable alternate option to the standard. Earth-like 
temperature restoration field. These modifications permit a 
tidally locked forcing scenario to be implemented, as is nec- 
essary for hot Jupiter type planets, with extreme tempera- 
ture differences between the star-facing dayside and shielded 
nightside. ( [Menou Rauscher|[2009 J . 



Teq(A,(/), Cr) 



where 



:Te7*(<T)+/?t.op(<j)AT(,(A, 



(2) 



- trop 



^stra + 



Zstva 
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Figure 1. Latitude- longitude plot of forcing temperature for the 
T42L15 run at cr = 0.7, showing the dayside 'hotspot' and cool 
nightside. The substellar point is at the centre of the image, and 
plots for other resolutions are identical. 



. ^ f sin (f (cr - Crstra)/(1 - CTstra)) CF > 



/5trop ( 

and 

AT0(A, 0) = cos(A) cos(0)ATep 



stra 
CTstra 



(3) 
(4) 

(5) 



The temperature produced by radiative-convective 
equilibrium, without winds or other factors, is represented 
by Teq, with the purely vertical part of this structure rep- 
resented by T^q^* and the purely horizontal part by /STq. 



The vertical element of the profile was chosen by|Menou 
Rauscher (2009[) to match that calculated by 



& Guillot 



Iro, Bezard 



2005) for HD 209458b. Such radiative-convective 
equilibrium temperatures may also be computed analytically 
from first principles using models such as that of | Guillot | 
(2010), which work has been extended to include the ef- 
fects of scattering by Heng et aL| ( |2012| . A scaling factor 
/3trop is applied to steadily decrease the temperature dif- 
ference between the dayside and nightside until it becomes 
zero above the tropopause, represented by crgtra and Zstra- 
Similarly, the temperature increment at the tropopause is 



given by ^Tstra (Menou & Rauscher 2009). The dry adia- 



batic lapse rate is represented by Ftrop, and the mean 'sur- 
face' temperature (temperature at the base of the model 
atmosphere) by Tsurf, with the equator-to-pole temperature 
difference denoted by ATep . Examples of the resulting tem- 
perature forcing pattern can be seen in Fig. [l] and Fig.|2] A 
'cold spot' is produced on the nightside of equal magnitude 
to the dayside hotspot. As a result, the zonal mean, or lon- 
gitudinally averaged, forcing temperature is the same from 
equator to pole on every level of the atmosphere, producing 
the vertical structure shown in Fig. [2] 

Runs were carried out at T42, T63, and T85 resolutions, 
with 15, 20, and 20 levels respectively, linearly spaced in cr. 
A full list of model, planetary, and run-specific parameters 
can be found in Tables [l] [2] and [3] respectively. All param- 



eters are taken from Menou & Rauscher (2009), who chose 
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Figure 2. Zonal mean plot of forcing temperature. The longitude 
averaging causes the day-night temperature difference to cancel 
out, so that the average equatorial temperature is equal to that 
at the poles. 



Table 1. Table of model parameters 



Parameter 


Symbol 


Value 


Ndel 




8 


Dissipation time-scale / day 


^diss 


5 X 10-3 


Rayleigh friction time-scale / day 


^frc 


oo 


Newtonian relaxation time /day 


^rest 


0.5 



HD209458b, a well-studied inflated hot Jupiter, and |Heng| 
■(|2011a|. 



et al. 



The hyper-diffusion time-scale on the smallest resolved 
scale is given by tdiss , while Ndel gives the order of that diffu- 
sion. As a result of the conservation of energy built into these 
models, energy 'builds up' at the resolution limit, requiring 
artificial diffusion to be applied to damp down and remove 
this effect, which would otherwise introduce unwanted ef- 
fects into the results. It is important to note that the hyper- 
diffusion time-scale applies to the smallest resolved scale in 
each run, and thus results in a weaker dissipation on any 
given scale for a higher-resolution run as compared with 
a lower-resolution one, effectively producing weaker overall 
diffusion. The effects of this can clearly be seen in Fig. |16| 
where increasing resolution permits the model to represent 
an energy cascade to smaller and smaller scales before being 
cut off by diffusion. The order and magnitude of the hyper- 
diffusion values are non-physical to an extent and require 
tuning to the simulation being run. A more detailed discus- 
sion of the effects of forcing and dissipation, with particular 
relevance to simulations such as that carried out here, can 



the planetary parameters to correspond to parameters for 



be found in "Thrastarson & Cho ( 2Qll| ). 

Rayleigh friction applies an effective frictional force to 
drag the winds towards zero velocity on each level on which 
it is set, using the time-scale tfrc; a time-scale of infinity 
thus indicates that no friction is applied. The Newtonian 
relaxation time, trest, determines the rate at which the model 
temperature is forced towards the input radiative-convective 
equilibrium state Teq. 

The 'surface', or base of the model, is defined as be- 
ing at the 1-bar pressure, Pq. The gravitational constant 
throughout the model is given by ^p, and the specific gas 
constant and heat capacity at constant pressure by Rs and 
Cp respectively. 
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Table 2. Table of planetary parameters 



Parameter 


Symbol 


Value 


Planetary radius / m 


a 




Rotation rate / 10~^ rad s"-*^ 




2.1 


Gravity / m s~^ 




8.0 


'Surface' pressure / bar 


Po 


1.0 


'Surface' temperature / K 


^surf 


1600 


Equator-pole temperature difference / K 


ATep 


300 


Tropopause temperature increment / K 




10 


Tropopause height / 10^ m 


^stra 


2 


Adiabatic lapse rate / 10~^ K m"-*^ 


Ptrop 


2 


Specific gas constant / J kg"-*^ K"-*^ 


Rs 


3779 


Heat capacity / J kg"-*^ K"-*^ 


Cp 


13226.5 


Table 3. Table of run-specific parameters 


Parameter Symbol T42L15 


T63L20 


T85L20 


No. of latitudes A^iat 64 


96 


128 


No. of longitudes A^ion 128 


192 


256 


No. of vertical levels A^iev 15 


20 


20 



Each run was carried out for a minimum of 350 plane- 
tary sidereal days following the determined spinup period, or 



approximately 1,000 Earth days (Heng et al. 2011a). Records 
were made ten times per day, equally spaced from one an- 
other, and all time averaging periods are over the 350 days 
unless otherwise stated. In this paper, the term 'day' shall 
always mean planetary sidereal days for consistency, noting 
that |Menou Rauscher|( |2009| ) and |Heng et"aLl ( |2011a) ) differ 
in their use of the term 'day', the former to mean planetary 
sidereal, the latter to mean Earth solar. 



3 RESULTS 

Fig.[3]shows the time mean temperature field in colour, with 
standard deviation contours overlaid. The model level clos- 
est to cr = 0.7 was chosen in each case, and all time means are 
taken over the 350-day period covering model days 30-380. 
It can be seen that the magnitude of both mean temperature 
and standard deviation is approximately consistent between 
runs, with standard deviation increasing slightly with res- 
olution, in particular noting that the T85L20 run in plot 
(c) shows no 15 K standard deviation contour towards the 
poles. The differences in the shape of the standard deviation 
contours between T42L15 and T63L20 may also be partially 
due to the different locations of the model levels, at a = 0.7 
and a = 0.725 respectively, an unavoidable result of using 
linear sigma spacing with different numbers of levels. The 
highest variance is found at ±30° N and approximately 100° 
E of the substellar point, where the equatorial jet has car- 
ried warm air from the dayside to the cool nightside, making 
the fluctuations caused by its north-south movement most 
apparent. Though small differences are apparent, all runs 
demonstrate very similar mean temperature fields and stan- 
dard deviations. The maximum mean temperature is 1698 
K at T42, 1705 K at T63, and 1709 K at T85, with corre- 
sponding maximum standard deviations of 67 K at T42, 64 
K at T63, and 66 K at T85. 
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Figure 3. Temperature mean and standard deviation for different 
resolution runs. Colours denote mean temperature, contour lines 
the standard deviation in Kelvin. The substellar point is at the 
centre of the image, or (0,0). (a) T42L15 at a = 0.7, (b) T63L20 
at cr = 0.725, (c) T85L20 at cr = 0.725. 



Versions of the T42L15 latitude-longitude instanta- 
neous 'snapshot' plots in fig. 3 of [Menou Rauscher (2009) 
are also shown for comparison in Fig. |4] While they do not 
precisely replicate the pattern of those images, the same gen- 
eral features can be seen, with the strong equatorial jet and 
weaker reverse flow beyond ±30° N producing the tempera- 
ture distribution seen in Fig. [4] While the wind s appear sim- 
ilar, the lack of a scale in Menou Rauscher) ( 2009 ) means 
that their strength cannot be directly compared in this im- 
age. Instead, reference must be taken from Fig. [S] which 
shows very similar zonal wind speeds to fig. 4 of [Menou fc| 
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Figure 4. 'Snapshot' of T42L15 temperature and wind vectors 
for comparison with [Menou &: Rauscher] ( |2009| >. The substellar 
point is at the centre of the image, (a) a = 0.7, (b) a = 0.37. 



Rauscher] ( |2009 |). Although the mean state is robust, the 
precise form of the equatorial jet and the temperature dis- 
tribution at a given moment is unpredictable, as the model 
is highly non-linear and thus its evolution is extremely sensi- 
tive, rendering such plots of limited use for the comparisons 
required in a benchmark test. 

The zonal mean temperature field, also averaged over 
the 350-day time period, is shown in Fig.|5] Note the differ- 
ent temperature scale to that of Fig. [3] due to the different 
temperature range encompassed. The most notable depar- 
ture from the forcing state's simple vertical profile can be 
seen in the centre of the figure, between ±40° N, where it 
can be seen that two 'hot spots' have developed, correspond- 
ing to the location of the temperature peaks seen either side 
of the equator in Fig. |3] These plots show the greatest dif- 
ference in standard deviation, with the location of the con- 
tours differing noticeably, particularly towards the top and 
bottom of the model atmosphere, although this may to some 
extent be due to the limited resolution in the vertical. All 
runs are noticeably cooler towards the poles below around 
a = 0.8 than those in fig. 6 of Heng et al. (2011a), which 
show temperatures of almost equal magnitude to those at 
the equator. 

By subtracting the time mean temperature field from 
the known forcing state and dividing through by the restora- 
tion time-scale, plots of the heating rate resulting from the 
forcing can also be derived, as shown in Fig. [6] It can be more 
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Figure 5. Time-averaged zonal mean temperature plots at differ- 
ent resolutions. Colours denote mean temperature, contour lines 
the standard deviation in Kelvin, (a) T42L15, (b) T63L20, (c) 
T85L20. 



clearly seen from this figure that average net cooling is thus 
experienced in the equatorially symmetric regions between 
0° and ±45° N, while net heating is experienced in regions 
poleward of ±45° N. Heating is seen in all locations below 
cr = 0.9, as well as between ±30° N above a — 0.2. The 
regions of negative heating correspond to the two warm re- 
gions noted previously in Fig. [5] The cooling regions appear 
to grow broader and the heating regions stronger, particu- 
larly towards the base of the model, with increasing resolu- 
tion. 

The plots in Fig. |7| display the time- aver aged, zon- 
ally averaged zonal (west-east) wind. Positive contours show 
wind directed west-to-east, 'out of the paper', negative con- 
tours the opposite. A strongly superrotating equatorial jet is 
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Figure 6. Time- averaged zonal mean heating rate plots at dif- 
ferent resolutions. Note the units of 10~^ K s-^ (a) T42L15, (b) 
T63L20, (c) T85L20. 



Figure 7. Time average of the zonal mean zonal wind. Positive 
contours and colours show wind directed 'out of the paper', neg- 
ative contours and colours the reverse, (a) T42L15, (b) T63L20, 
(c) T85L20. 



seen between 25° S and 25° N and below a — 0.2, with peak 
windspeeds of approximately 1,200 m in all cases, and 
there is a corresponding broader and weaker return flow out- 
side this region. Maximum and minimum mean windspeeds 
are found to be 1220, -687 m s"^ for the T42L15 run, 1200, 
-693 m s"^ at T63L20, and 1179, -698 m s"^ at T85L20, 



comparable to the Heng et al. (2011a) 'shallow hot Jupiter' 



spectral model using the same parameters. From these re- 
sults, differences of order 1% are noted when varying the 



resolution and hyper diffusion. |Heng et aL (2011a) note that 
mean winds are uncertain at the level of approximately 10% 
in their 'deep hot Jupiter' simulation when changing pa- 
rameters. The jet is unbounded at the frictionless base of 
the model, but closes towards the top at around a = 0.2, 



with weak reverse flow above. This pattern corresponds well 
to that seen in both Menou Rauscher] ( 2009 ) and the 



Heng et al. ( 2011a| ) paper for the 'shallow hot Jupiter' test 



case, although again, Menou & Rauscher (2009 ) display only 



snapshots. The corresponding zonal wind snapshot for the 
T42L15 run is shown in Fig.|8] and is broadly similar to that 
of Menou Rauscher^ (2009) fig. 4, with the most notable 



difference being that the greatest zonally- averaged westward 
windspeed is approximately 830 as opposed to 730 m s~^ at 
this particular instant in time. 

The time evolution of the temperatures and wind speeds 
at various representative areas of the T42L15 run up to day 



100, for comparison to fig. 5 of Menou & Rauscher (2009 ), is 



© 0000 RAS, MNRAS 000, 000-000 



Benchmark experiments for exoplanetary GCMs 7 




0° 

Latitude 



-1200 -800 



-400 400 
Zonal wind / m s ' 



800 1200 



Figure 8. T42L15 zonal wind 'snapshot' plot, for comparison 
to fig. 4 of Menou & Rauscher (2009). Positive colours and con- 
tours show wind directed 'out of the paper', negative colours and 
contours the reverse. 




CD 

E 1450 

(D 

1400 
1350 



40 60 
Time / days 



100 



Figure 10. Representative T plots for the T42L15 run, as in fig. 
5 of [Menou Rauscher] ( [2009| >. Light and dark red lines repre- 
sent temperatures at the sub- and anti-stellar points, respectively, 
green and dark green dotted lines temperatures at the east and 
west limbs respectively, and blue and dark blue dashed lines tem- 
peratures at the north and south poles, respectively. 
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Figure 9. Representative U plots for the T42L15 run, as in fig. 
5 of [Menou Rauscher] ( pOOQl . Solid line represents zonal mean 
wind at the equator, dashed line the equatorial maximum. 



shown in Fig. |9] and Fig. [To] Rather than use single points, 
representative areas were chosen, averaging around the lon- 
gitude circles at ±87.9° N for polar values, and over the 
latitude band between ±4.2° N for equatorial values, which 
were also averaged over the longitude points between ±2.8° 
of the specified location (e.g. the substellar point). This en- 
ables greater reproducibility, and also helps to avoid model- 
specific issues in choosing the points: PUMA, for example, 
has neither points at ±90° N nor exactly 0° N. For all runs, 
the data were sampled ten times per day. The plots shown 
demonstrated extremely high variability over small time- 
scales, and were smoothed using a simple boxcar smooth of 
width 1 day, or 11 records, to better display overall trends. 
Such high variability shows a requirement for high-frequency 
sampling in the creation of these plots, as daily sampling 
does not capture the full extent of the variability and can 
miss high-frequency features altogether, or produce spurious 
signals through aliasing. 

Fig. 11 (a), (c), and (e) display the streamfunction on 
the model level closest to a = 0.7. This demonstrates the cir- 



culation of the flow, and its direction, with a large, positive 
streamfunction value indicating strong clockwise circulation. 
Its associated variance can be seen in Fig. 11 (b), (d), and 
(f). It can be seen that there are two major circulation fea- 
tures, one to either side of the equator, centred around 60° E 



of the substellar point, while the greatest variability is found 
on the equator around 60° W of the substellar point. The 
circulation at this level strengthens slightly with increasing 
resolution, with its maximum value just exceeding 4 x 10^° 
m^ s"^ at T85. 

Fig. [12] shows snapshots of vorticity, a measure of the 
local rotation at each point, on the level nearest a = 0.7 for 
both T42L15 and T85L20 runs, demonstrating the impact of 
increasing resolution. Although a degree of structure is visi- 
ble at T42, much finer, smaller-scale structure can be seen in 
the T85 plot, with long 'streamers' of high-magnitude vor- 
ticity visible that are washed out at T42 resolution. A single 
time mean is shown as the time mean vorticity plots are 
almost identical between resolutions. Though fields such as 
temperature and wind appear quite similar to one another 
even in snapshot forms, the vorticity clearly shows the ben- 
efit of higher-resolution runs, enabling much finer detail to 
be resolved. 

The plot in Fig. [13] displays the mean meridional circu- 
lation, a measure of the mass of air circulating about a given 
point. Positive contours indicate clockwise circulation, neg- 
ative anticlockwise; a typical MMC plot for the Earth would 
show positive contours between 0°N and 30° N, illustrat- 
ing clockwise circulation in which air rises over the equator 
and descends further north, followed by the inverse between 
30° N and 60° N, and a further clockwise circulation between 
there and the pole ([ Peixoto Oort|1992 ) . Similarly, the in- 
verse pattern is seen in the southern hemisphere. Here it can 
be seen that there are two main circulation features, with 
air descending rather than rising over the equator (as ex- 
pected from Showman Polvani|20lT ), and rising between 
±30° N and ±60° N, with a weak reverse circulation towards 
the poles. The equatorial circulation contracts and weak- 
ens slightly with increasing resolution, with the two distinct 
peaks at different altitudes becoming more obvious 

Fig. [M] shows the local superrotation index, which is a 
measure of the degree to which the angular momentum of 
each element of the atmosphere exceeds that which it would 
have in solid-body rotation. The local superrotation index s 
is defined by 



(6) 
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Figure 11. Streamfunction mean (left) and standard deviation (right) for different resolution runs. Plots (a) and (b): T42L15 
at cr = 0.7; (c) and (d): T63L20 at cr = 0.725; (e) and (f): T85L20 at cr = 0.725. The substellar point is at the centre of the 
image. Mean contours: lO-'^^m^s"-'^ ; standard deviation contours: lO^m-^s"-*^ 



where m is the axial angular momentum per unit mass of 
the atmosphere derived from the zonal mean zonal wind 
the longitudinal average of the zonal (east- west) component 
u of the wind field (Lewis & Read 2003). In general, the 



found inlShowman & Polvanil (12011 1). Only the T85L20 run 



axial angular momentum per unit mass is given by 



m — Qa^ cos^((/)) + ua cos{(j)) 



(7) 

A global superrotation index S can also be calculated by 
integrating over the whole atmosphere: 



(8) 



[ma cos{(j)) / g)d\d(j)dp \ /Mo — 1 



where Mo is the same integral for the atmosphere at rest, or 
1^ = 0. 

A westerly wind (blowing west-to-east) over the equa- 
tor cannot be created from an atmosphere initially at rest 
simply by moving air parcels from other regions of the at- 
mosphere, since the maximum angular momentum available 
is that at the equator. The existence of superrotation is thus 
a signature of eddy processes occurring in the atmosphere, 
transporting angular momentum equatorward. A detailed 
study of superrotation under hot Jupiter conditions can be 



is displayed, as the results for each run are visually identical. 
The maximum value lies between 0.56 and 0.57 in each case, 
while the minimum is —1 at the poles. 

Detailed study of the angular momentum budget over 
the course of each run reveals that the global superrotation 
index begins at 5* = 0, as expected from the model's initial- 
isation state of zero wind. It then climbs over the first 15 
days to a value of 0.033 =b 0.003 in each simulation, indicat- 
ing that angular momentum is not fully conserved, and ad- 
ditional energy has been imparted to the atmosphere. With 
no diurnal tides, surface friction, or topographical features 
to provide this extra momentum, it is likely to have been 
acquired through model dissipation. 

Fig. [15] shows the time evolution of the global statis- 
tics RMS vorticity and RMS divergence, which are output 
at each model timestep, over the initial 400-day period of 
the T42L15 run. These values are the root mean square of 
the entire global vorticity and divergence fields. A stable 
state is reached after 25-30 days, after which the observed 
degree of variability is unchanged; the same is true at the 
other resolutions, although the mean values reached differ. 
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Figure 12. Vorticity day-100 snapshots and time mean for T42L15 and T85L20 runs, (a): T42L15 snapshot at cr = 0.7; (b): 
T85L20 snapshot at cr = 0.725; (c): T85L20 time mean at cr = 0.725. The substellar point is at the centre of the image. Note 
the difference in scale of plot (c), at 10~^ s"-*^ rather than the lO"'^ s"-*^ of plots (a) and (b) 



As with the representative area plots of Fig. [9] and Fig. |1Q| 
this demonstrates that the spin-up time for puma under 
these conditions is approximately 25-30 days, after which 
data may be taken without fear of compromising the results 
with spurious spin-up effects. 

Fig. [16] shows the kinetic energy spectra for three 
resolutions. Notably, low (large-scale), even (symmetric) 
wavenumbers have much higher amplitudes than their odd- 
valued counterparts, due to the highly equatorially sym- 
metric, large-scale nature of the thermal forcing and final 
state. The dotted line has a slope of —3, that expected 
from an enstrophy-cascading range in two-dimensional tur- 
bulence ( |Kraichnan| [l967 ) . The majority of the spectrum 
lies closely parallel to this line, demonstrating that this 
regime holds over most modelled scales. This quasi- two- 
dimensional regime is to be expected from the scales reach- 
able by these studies, as the smallest resolved scale even 
at maximum resolution, T85, is still on the order of 10^ 
km, with flow on this scale strongly constrained by the ef- 
fects of planetary rotation and atmospheric depth, rather 
than fully three-dimaensional turbulence ( HoughtQii||1986 ) . 
Higher wavenumbers correspond to smaller scale features, 
and the greater kinetic energy present at higher wavenum- 
bers in the higher resolution runs thus results in the greater 
detail and higher extrema seen most clearly in the vorticity 
plots of Fig. [12] In each case, the spectrum tails off sharply 
towards the run's wavenumber cutoff, with a slope of around 
— 15. This sharp decrease near the cutoff is not linked to 
physical expectations and is a result of the model diffusion. 



While diffusive processes do naturally occur, the limitations 
of modelling require that they must be represented at pro- 
gressively larger scales (lower wavenumbers) as the model 
resolution decreases, to avoid an unphysical build-up of en- 
ergy at the smallest resolved scales. In a true system, this 
energy would continue to cascade down to ever smaller scales 
and eventually be dissipated; the model system, however, 
cannot reach such small scales, and since it conserves en- 
ergy efficiently, must have additional dissipation applied. 



4 DISCUSSION 

As the results presented have shown, although 'snapshot' 
plots are useful for capturing a view of how the simulated 
atmosphere is behaving, they are effectively irreproducible 
(contrast, for example. Fig. |4] of this paper and fig. 3 of 
Menou Rauscher (2009)), and can only be compared qual- 
itatively, rather than quantitatively. They may even be mis- 
leading, if the 'snapshot' time is chosen while the model is far 
from the mean state. Time mean plots, as produced by |Heng| 
et al. (2011a), provide directly comparable results, as the 



mean flow would be expected to be very similar between any 
two simulations reaching the same stable state. However, the 
variability of the model is then lost without the addition of 
variance or standard deviation information. Such small-scale 
variability is another important point of comparison, as it 
demonstrates potential planetary 'weather' and also shows 
the mean transport of such quantities as heat and momen- 
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Figure 13. Mean meridional circulation. Positive contours cor- 
respond to circulation in a clockwise sense, negative contours to 
anticlockwise circulation, (a) T42L15, (b) T63L20, (c) T85L20 



turn by transient waves. Taken together, therefore, plots of 
the mean and associated standard deviation are found to be 
more useful for the purposes of model intercomparison work. 

Resolution is found to make little difference to many of 
the mean results studied, with similar mean temperatures 
and winds recorded in all cases. The most visible differences 
between runs are seen in the standard deviation contours, 
denoting different levels of variability in different regions. 
The persistence of the large-scale features between runs and 



models (Menou & Rauscher 2009 Heng et al. 2011a) sug- 



gests that such features, with a powerful superrotating equa- 
torial jet and correspondingly offset temperature maxima 
and minima, are likely to be observed on true extrasolar 
planets subject to such extreme forcing conditions. The ob- 
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Figure 14. Time-averaged local superrotation index for the 
T85L20 run. 
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Figure 16. Kinetic energy spectra for all runs. The spectra have 
been averaged over all levels and over a period of ten days. The 
blue lines denote the spectrum of the T42L15 run, the red lines 
T63L20, and the green lines T85L20. The total energy spectrum 
is shown as a solid line, while the zonal component is dashed and 
the eddy dash-dotted. The dotted line is a reference line with a 
slope of —3. 
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servations available appear to support this conclusion, ev- 



idenced by the results of Knutson et al. (2007), who fit a 



simple model of planetary brightness temperature to the 
phase variation of the HD 189733 light curve to obtain a 
temperature map of HD 189733b, [Majeau, Agol &: Cowan| 
( |2012| , who use two different methods to gain a temperature 
map of HD 189733b from its secondary eclipse, and Snellen 
et al. (2010), who detected a 2 km s~^ CO blueshift imply- 



ing strong winds flowing from the dayside to the nightside 
of HD 209458b. However, very few such observations have 
been made, and sufficient resolution to confirm or deny the 
presence of features at scales much below the global is still 
unattained. For the present, model results thus remain the 
only available avenue of study for the majority of features. 

The analysis of Fig. [15] as well as Fig. [5] and Fig. [lO] 
shows that the approximate spin-up time for puma is around 
25-30 days, after which a stable state is reached and the 
starting conditions are effectively erased. This determines 
the number of records it is necessary to discard before av- 
eraging the data to avoid introducing spurious information, 
and may vary from model to model. This is a relatively short 
period of time compared to typical Earth spinup times, and 
extremely short in comparison to the times required to spin 
up a model of Jupiter itself, which receives very little exter- 
nal heating. 

Sampling frequency can visibly alter the appearance of 
the representative value plots, with the daily sampling used 



by Menou & Rauscher ( 2009 ) hiding the high variability seen 



on smaller time-scales, and use of a time filter together with 
a high sampling frequency is recommended. 

The kinetic energy spectra at all resolutions largely fol- 
low the —3 law of Kraichnan (1967), but demonstrate a 



sharp fall in energy at scales approaching the truncation 
wavenumber, with a slope of —15 or steeper. Each increase 
in resolution studied produces an extension of the enstrophy- 
cascading range before the point at which dissipation takes 
over is reached. Thrast arson & Cho ( [2011 ') suggest that 
model calculations with these large forcing amplitudes and 
short restoration and dissipation time-scales are likely to be- 
come over- or under-damped, and the use of energy spectra 
and fields such as vorticity are recommended to determine 
whether this is the case. 



5 CONCLUSIONS 

Benchmark tests of GCMs for hot Jupiters have been con- 
sidered, and a variety of diagnostics produced and analysed. 

'Snapshot' plots are of limited use for the purposes of 
an intercomparison study because the precise phase of waves 
depends sensitively on the initial conditions and the evolu- 
tion of the flow. The use of mean and standard deviation 
plots is therefore encouraged. The addition of plots such as 
streamfunction, mean meridional circulation, and superro- 
tation index to the simple temperature and wind fields are 
also suggested, as these alone do not capture all aspects of 
the simulation, although more care is required in interpret- 
ing some such plots. 

The use of high sampling frequencies is also found to 
be important, with the degree of variability seen in the 'rep- 
resentative area' plots of Fig. [5] and Fig. [lO] demonstrably 
dependent on the frequency at which the data is sampled. 



A sampling frequency of at least ten records per day is sug- 
gested, and the data time-filtered. 

The resolution of the model runs is found to have little 
effect on the overall mean and variability of the data. How- 
ever, much finer detail is captured at high resolution, as illus- 
trated in the vorticity plots of Fig. [12] Smoother fields such 
as temperature show less visible effect of this increased res- 
olution. The energy spectra show that the higher-resolution 
models continue to follow the —3 law to smaller scales be- 
fore dissipation takes effect, demonstrating that worthwhile 
data can be gained by increasing the resolution. 

The large-scale features of the temperature and wind 
fields correspond in form and magnitude to those of |Menou] 
Rauscher (2009) and Heng et al. (2011a), to the extent 



that those can be determined from the results provided. This 
reinforces the growing consensus that such features, with a 
'hotspot' offset from the substellar point by a strongly su- 
perrotating equatorial jet, are likely to be found on a variety 
of hot Jupiter type exoplanets. 
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